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Abstract. We study the dynamics of water confined between hydrophobic flat 
surfaces at low temperature. At different pressures, we observe different behaviors 
that we understand in terms of the hydrogen bonds dynamics. At high pressure, the 
formation of the open structure of the hydrogen bond network is inhibited and the 
surfaces can be rapidly dehydrated by decreasing the temperature. At lower pressure 
the rapid ordering of the hydrogen bonds generates heterogeneities that are responsible 
for strong non-exponential behavior of the correlation function, but with no strong 
increase of the correlation time. At very low pressures, the gradual formation of the 
hydrogen bond network is responsible for the large increase of the correlation time and, 
eventually, the dynamical arrest of the system and of the dehydration process. 
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1. Introduction 

1.1. Supercooled Water and Polyamorphism 

Supercooled water is a metastable state of the liquid water with respect to the crystal 
ice and occurs when H2O is cooled slowly below the freezing point. It can be easily 
prepared by cooling relatively pure commercial water and is commonly observed in 
nature, especially under confinement as for example in plants where water is still liquid 
at — 47°C. It may crystallize into ice by catalysts, such as an ice crystal, or by a large 
exchange of mechanical energy, such as a strong vibration. However, without catalysts 
or external work, it is a useful example of a liquid whose dynamics slows down when 
the temperature decreases toward the homogeneous nucleation temperature (with 
Th = -41°C at 1 atm and Th = -92°C at 2000 atm) p. The study of supercooled 
water is important in the long-lasting efforts to understand the physics of the hydrogen 
bonds, which characterize water and are the origin of its unique properties [2]. 

Below Th water crystallizes, but a fast quench from ambient T toT < Tx = — 123°C 
(at 1 atm) freezes the liquid into a disordered (not crystalline) ice, leading to the 
formation of (amorphous) glassy water [3]. By increasing the pressure it has been 
observed that glassy water changes from low density amorphous (LDA) to high density 
amorphous (HDA) [4] and from HDA to very-high density amorphous (VHDA) [5]. 
Therefore, water shows polyamorphism and this phenomenon has been related to the 
possible existence of more than one liquid phase of water [6]. Other thermodynamically 
consistent scenarios have been proposed, such as a first-order liquid-liquid phase 
transition with no critical point [7] or the absence of liquid-liquid phase transition itself 
[H], and their experimental implications are currently debated 

1.2. Confined Water 

Confinement is an effective way to suppress the crystal nucleation. Water molecules in 
plants, rocks, carbon nanotubes, living cells or on the surface of a protein, are confined 
to the nanometer scale. Many recent experiments show a rich phenomenology which is 
the object of an intense debate (see for example [IOl[IIl[l2l[T3l[71[Tl[T5l[l6l[I7l[l8l[T9]). 
The relation between the properties of confined and bulk supercooled water is not fully 
understood, therefore it is useful to elaborate handleable models that could help in 
developing a consistent theory. 

1.3. Models for Water 

Since the experiments in the supercooled region are difficult to perform, numerical 
simulations are playing an important role in the interpretation of the data. These 
simulations are based on models with different degree of detailed description and 
complexity. Among the most commonly used, there are the ST2 |20j (two positive and 
two negative charges in tetrahedral position plus a central Lennard- Jones or LJ), the 
SPC and its modification SPC/E (three negative and one positive charges plus a central 
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LJ) [2T], the TIP3P, TIP4P, and TIP5P with three, four and five charges, respectively 
[22] . the models with internal degrees of freedom [23], polarizability and three body 
interactions [23] , and ab initio models for the evaluation of the significance of quantum 
effects [25] . 

Models with increasing details lead to improved accuracy in the simulations, but 
suffer from important technical limitations in time and size of the computation. Due 
to their complexity, an analysis of how each of their internal parameter affects the 
properties of the model is in general very computationally expensive and a theoretical 
prediction of their behavior is unfeasible, limiting their capacity to offer an insight into 
unifying concepts in terms of universal principles and hampering their ability to study 
thermodynamics properties. 

On the other hand, simplified theoretical models, such as the Ising model, can be 
very drastic in their approximations, but are able to offer a general view on global, 
and possibly universal, issues such as thermodynamics and dynamics. They can be 
studied analytically or by simulations, providing theoretical predictions that can be 
tested by numerical computation. Therefore, they offer a complementary approach 
to the direct simulation of detailed models. Their characterization in terms of few 
internal parameters make them more flexible and versatile, offering unifying views of 
different cases. Their study by theoretical methods, such as mean fleld, leads to analytic 
functional relations for their properties, allowing to predict their behavior over all the 
possible values of the external parameters. Their reduced computational cost permits 
direct test of the theoretical forecasts in any condition. Furthermore, their study helps in 
deflning improved detailed models for greater computational efficiency. Therefore, there 
is a strong need for simplifled models of water that can keep the physical information 
and provide a more efficient simulation tool. 

In the following we will describe a simple model for water in two dimensions (2D). 
The model can be easily extended to 3D, however we consider here only the 2D case. It 
can be thought as a monolayer of water between two fiat hydrophobic plates where the 
water- wall interaction is purely repulsive (steric hard-core exclusion). It has been shown 
that the properties of confined water are only weakly dependent on the details of the 
confining potential between smooth walls and that the plate separation is the relevant 
parameter determining whether or not a crystal forms [26]. The case we consider here 
corresponds to the physical situation in which the distance between the plates inhibits 
the crystallization. 

2. Cooperative Cell Model Confined Between Partially Hydrated 
Hydrophobic Plates 

2.1. General Definition of the Model 

Cell models with variable volume are particularly suitable to describe fiuids at constant 
pressure P, temperature T and number of molecules as in the majority of experiments 
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with water. Here we consider a cell model based on the idea of including a three-body 
interaction, as well as an isotropic and a directional term for the hydrogen bonds [27|l28]. 
The three-body term takes into account the cooperative nature of the hydrogen bond. 
This model reproduces the experimental phase diagram of liquid water [29l [30] and, by 
varying its parameters, the different scenarios for the anomalies of water [311 [32] . 

We consider water confined between two hydrophobic plates within a volume V 
that accommodates N molecules. To minimize the boundary effects, we consider double 
periodic conditions. We assume that the distance between the plates is such that water 
forms a monolayer on a square surface, similar to the one observed in molecular dynamics 
simulations for TIP5P water, where each water molecule has four nearest neighbors |33j . 
The system is divided into cells of equal volume v = V/L'^. Differently from previous 
works on this model [34l[27l[35l[29l[2Hl[9l[2l[30l[3ll[32], here we consider partially 
hydrated walls with a hydration level corresponding to a water concentration of 75% 
with respect to the available surface of one wall. This concentration is well above the 
site percolation threshold (59.3%) on a square lattice [36|. We assume that in each cell 
i = 1, 2, ... there is at most one water molecule and we assign the variable rij = 1 if 
the cell is occupied, and Ui = otherwise. 

The water-water interaction is described by three terms. One takes into account all 
the isotropic contributions (short-range repulsion of the electron clouds, weak attractive 
van der Waals or London dispersion interaction, strong attractive isotropic component of 
the hydrogen bond interaction). This term is given by a hard-core repulsive at distance 

1/2 

volume Ro = Vq plus a LJ potential 



where r is the distance between two molecules. 

The second term takes into account the directionality of the hydrogen bonds 
interaction. Each molecule participates in at most four hydrogen bonds when each 
donor 0-H forms an angle with the acceptor O that deviates from linear (180°) within 
30° or 40° approximately. To incorporate this feature to the model, we associate to 
each water molecule i four indices aij G [1, 2, ... g], one for each possible hydrogen bond 
with the nearest neighbors molecules j. The parameter q is chosen by selecting 30° as 
maximum deviation from the linear bond, hence q = 180/30 = 6. With this choice, 
every molecule has q'^ = 6'^ = 1296 possible orientations. This directional interaction 
can be accounted for by the energy term 



where the sum runs over all nearest-neighbor cells, J > is the energy associated to 
the formation of a hydrogen bond, and Nhb is the total number of hydrogen bonds in 
the system. 




for r < R[ 



for r > Ri 
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As a consequence of their directionality, the hydrogen bonds induce a rearrangement 
of the molecules in an open network (tetrahedral in 3D), in contrast with the more 
dense arrangement that can be observed at higher T or higher P jST]. This effect can 
be included in the model by associating a small volume contribution vhbi due to the 
formation of each bond, to the total volume of the system 

V = Vo + NhbVhb (3) 

where Vq = Nvq is the volume of the liquid with no hydrogen bonds. 

Finally, experiments show that the relative orientations of the hydrogen bonds 
of a water molecule are correlated [38]. As an effect of the cooperativity between 
the hydrogen bonds, the 0-0-0 angle distribution becomes sharper at lower T [llj . 
Therefore, the accessible q values are reduced. We include this effect in the model by 
the energy term 

Ea = -Ja^UiY^ (4) 

where is the energy gain due to the cooperative ordering of the hydrogen bonds, 
{k, l)i stands for the 6 different pairs of bond indices (arms) of a molecule i. At low T 
this terms induces the reduction of accessible values of q (orientational ordering). In the 
limiting case J^j = there is no correlation between bonds [8]. This case corresponds 
to consider as negligible the observed change in 0-0-0 angle distribution with T. The 
case Jo- 1 corresponds to a drastic reduction of orientational states for the molecule 
from to q, with all the relative orientation between the hydrogen bonds kept fixed 
independently from T and P. 



2.2. The Parameters Dependence of Phase Diagram of the Model 

The complete enthalpy for the fluid reads 

H = U{r) + Ej + E^ + PV. (5) 

Each of the energy terms of ([5]) provides an energy scale, and depending on their 
relative importance different physical situations emerge. For instance, the choice J < e 
guarantees that the hydrogen bonds are mainly formed in the liquid phase. As showed 
in [32] , when = the model gives rise to the scenario with no liquid-liquid coexistence 
[8]; at larger Ja/J the model displays the scenario with a liquid-liquid critical point at 
positive pressure [6]; further increase of Ja/J brings to a liquid-liquid critical point at 
negative pressure |39]; finally more increase of J^/J brings to the scenario where the 
liquid-liquid phase transition extends to the negative pressure of the limit of stability 
of the liquid phase [7] . 
The term 

Ej + PV = -{J - Pvhb)Nhb + PVo (6) 

is proportional to the number of hydrogen bonds Nhb and the associated gain in 
enthalpy decreases linearly with increasing pressure, being maximum at the lowest 
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accessible P for the liquid, and being zero at P = J/vhb- For P > J/vhb the formation 
of hydrogen bonds increases the total enthalpy. 

In the following we choose J/e = 0.5, J^/e = 0.05, vhb/vo = 0.5 as in Ref. [29], 
for sake of comparison. In the next sections we describe the simulation method adopted 
here and the results. 

3. Kawasaki Monte Carlo Method 

Monte Carlo (MC) simulations are performed in the ensemble with constant A^, P, T 
for = 1875 water molecules distributed in a mesh of = 2500 cells, corresponding to 
75% water concentration. We do not observe appreciable size effects at constant water 
concentration for (A^ = 300, = 400) and {N = 1200, = 1600). Preliminary results 
at 90% water concentration display no significant differences with the present case. 

Since we are interested in studying the dynamic behavior of the system when 
diffusion is allowed, we adopt a MC method a la Kawasaki [IQ]. This method consists 
in exchanging the content of two stochastically selected nearest neighbors (n. n.) cells 
and accepting the exchange with probability 



where AH is the difference between the enthalpy of the final state and that of the initial 
state. 

The same probability is adopted for every attempt to change the state of any of 
the bond indices aij. One MC step consists of 5A^ trials of exchanging n. n. cells i and 
j or changing the state of the bond indices aij (every time selecting at random which 
kind of move and which i and j are considered) followed by a volume change attempt. 
In the latter case, the new volume is selected at random in the interval [V — SV, V + SV] 
with SV = O.bvo- 

For sake of simplicity, in the Lennard- Jones energy calculation we do not consider 
the volume changes due to the hydrogen bonds formation and we allow the formation 
of hydrogen bonds independently from the size of the n. n. cells. We include a cut-off 
for the Lennard- Jones interaction at a distance equal to three times the size of a cell. 
Despite with this choice the cut-off is density dependent, we tested that in this way we 
largely reduce the computational cost and that the results do not change in a significant 
way with respect to the case without the cut-off. 

In the simulations we follow an annealing protocol, by preparing a random 
configuration and equilibrating it at high T and low P (gas phase). We equilibrate 
the system over 5 x 10^ MC steps, and calculate the observable averaging over at least 
2 X 10^ MC steps. We use the final equilibrium configuration at a given temperature as 
the starting configuration for the next lower temperature. 
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4. Results and Discussion 

4-1. The Phase Diagram 

From the calculation of the density of the system at fixed P as function of T we find 
a liquid-gas first order phase transition ending in a critical point C at approximately 
kBTc/e = 1.9 ± 0.1, PcVo/e = 0.50 ± 0.05 (Fig. [1^). In the hquid phase, we observe 
a temperature of maximum density (TMD) for any isobar with Pvo/e < 1. This 
anomalous behavior in density is one of the characteristic features of water. Here we 
note that Pv^/e < 1 corresponds to the limiting value P = J/vhb above which the 
formation of bonds increases the enthalpy in Eq.(l6]), consistent with the intuitive idea 
that the density maximum in water is a consequence of the open structure due to the 
hydrogen bond formation. For P > J/vhb the formation of hydrogen bonds is no 
longer convenient in terms of free energy and the system recovers the behavior of a 
normal liquid (with no open structures). 

At lower T and for P < J/vhb = le/'i^o, the isobars display another discontinuity 
between high density and low density (Fig{T]D), both at p higher than the gas density. 
As in the case with no diffusion [29], we associate this discontinuity in p to a first order 
phase transition between a high density liquid (HDL) and a low density liquid (LDL) 
ending in critical point C at kBTc/e = 0.06 ± 0.01 and Pc'Vo/e = 0.60 ± 0.15. The 
values of the critical temperature and pressure for C and the critical temperature for C 
are consistent, within the error bars, to those previously estimated [29l [28] , while the 
value of Pc' is somehow smaller. This is due to the difficulty in estimating the position 
of C from the discontinuity of p and it will be discussed in a future work. Here, instead, 
we focus on the analysis of the dynamics of the system at low T. 

4-2. The Dynamics at Low Temperature 

To get an insight into the processes at low T associated with the hydrogen bond 
dynamics, we calculate the the quantity Mi = ^'^jO'ij, that gives the orientational 
order of the four bond indices of molecule i, and its correlation function 



where the (■) stands for the thermodynamic average [30l [3T]. 

As we have discussed in the previous section, the effect of the hydrogen bond is 
evident only at pressure P < J/vhb and below the TMD line. Therefore we focus on 
state points at P < le/vQ and T < OAe/kB- 

• We first observe that for P = le/vQ the correlation function is always exponential 
for all the temperatures with 0.005 < kBT/e < 0.4 (Triangles in Figj2]). From 
the exponential fits (dotted lines in Fig 12]) the characteristic decay time is r ~ 1.3 
MC Steps (MCS) at kBT/e = 0.4 (green symbols), r ~ 2.2 MCS at kBT/e = 0.1 
(red symbols), r ~ 10 MCS at kBT/e = 0.05 (black symbols). Therefore, the 
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correlation time increases of one order of magnitude for a decrease of T of one 
order of magnitude, but r is always much shorter than the simulation time. 
This behavior is consistent with the observation of the typical configurations of 
hydrogen bonds at these P and T (upper row in Fig. [3]). The water molecules 
condense on the hydrophobic surface, leaving a large dehydrated area, but the 
number of hydrogen bonds does not increase in an evident way and no orientational 
order is observed even at the lowest T. 

Next we decrease the pressure to P = 0.7e/vo (Diamonds in Fig|2]), close to the 
the liquid-liquid critical pressure. At ksT/e = 0.4 (green symbols), the correlation 
function is exponential with r ~ 1 MCS. At lower T the correlation function is no 
longer exponential. It can be well described with a stretched exponential function 

CMit)=Coexp[-it/rf'j (9) 

where Cq, t and (3 < 1 are fitting constant. For (3 = 1 the function is exponential, 
and the more stretched is the function, the smaller is the power < (3 < 1. 
At kBT/e = 0.1 (red symbols), the stretched exponential (continuous red line) has 
f3 = 0.8 and r ~ 3.1 MCS, while at ksT/e = 0.05 (black symbols), the stretched 
exponential (continuous black line) has j3 = OA and r ~ 5 MCS. 
Hence, along isobars close to the liquid-liquid critical pressure, the dynamics of the 
liquid is not much slower than that at Pvo/e = 1, but it becomes non exponential 
below the TMD, a typical precursor phenomenon of glassy systems |Hj. This 
phenomenology is usually associated with the presence of heterogeneity in the 
system. While in a homogeneous system the characteristic relaxation time is 
determined by some characteristic energy scale, in a heterogeneous system the 
presence of many characteristic energies induces the occurrence of many typical 
time scales, giving rise to a non exponential behavior. 

This is consistent with the increasing number of hydrogen bonds that can be 
observed at this pressure in the high density liquid phase (central panel on the 
second row in Fig. [3]). The formation of many hydrogen bonds inhibits the free 
diffusion of the molecules, determining heterogeneities. 

This effect is extremely evident at the lower ksT/e = 0.05, where the non 
exponential behavior of the correlation function is very strong {(3 = 0.4) and the 
hydrogen bond network is fully developed, but with many isolated locally disordered 
regions (left panel on the second row in Fig. [3]). This is the typical configuration 
of the hydrophobic surface at low T at this P. The surface is partially hydrated 
by low density liquid water forming an open network of bonded molecules. The 
network is almost complete, but many small heterogeneities are present. 

By further decreasing the pressure to P = 0.4e/t>o (Squares in Figl2]), we observe 
an exponential decay of the correlation function at ksT/e = 0.4 (dotted line) with 
r ~ 1.3 MCS and a stretched exponential behavior when the liquid-liquid critical 
temperature is approached. Differently from the case at P = 0.7e/vo, at P = OAe/vo 
the stretched exponent does not decrease in a strong way, but the relaxation time 
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largely increases at lower T. 

In particular, at ksT/e = 0.1 (red symbols), the stretched exponential (dashed red 
line) has (3 = 0.9 and r ~ 18 MCS, while at ksT/e = 0.05 (black symbols), the 
stretched exponential (dashed black line) has (3 = 0.7 and r ~ 100 MCS. Therefore, 
the correlation time slows down of two order of magnitudes by decreasing the T 
from keT/e = 0.4 to keT/e = 0.05. 

This different behavior with respect to the case at P = 0.7e/vo could be associated 
to the different evolution of the hydrogen bond network. Previous mean field [28] 
and MC results [2] have shown that the increase of the number of hydrogen bonds 
with decreasing T is much more gradual at low P than at high P. Therefore, at 
high P this process could generate quenched-in heterogeneities and strong non- 
exponentiality. At low P, instead, the gradual building up of the hydrogen bond 
network freezes the system, inducing a large slowing down. As a consequence, the 
systems has large relaxation times, but the hydrogen bond network develops in a 
more ordered way (third row in Fig. [3]). 

The gradual increase of the number of hydrogen bonds is also responsible for the 
slower dehydration of the hydrophobic surface. Many small dehydrated regions 
seem to be annealed at a very slow rate by the progressive network formation. 

• We finally consider the lowest pressure P = O.le/vo (Circles in Figj2]). In this 
case, the correlation time is exponential for kBT/e = 0.4 (dotted green line) and 
ksT/e = 0.1 (dotted red line), with r ~ 1.3 MCS and r ~ 2.2 MCS, respectively. 
At these temperatures the system is homogeneous, with no order at ksT/e = 0.4 
(right panel on the bottom row in Fig. [3l) and with small orientationally ordered 
regions at ksT/e = 0.1 (central panel on the bottom row in Fig. [3]). These regions 
are the effect of the slowly increasing number of hydrogen bonds Nhb and are 
possible only at low P, because only at low values of P the enthalpy increase 
proportional to Nhb [Eq-dS])] is small. 

As a consequence of the gradual increase of these locally ordered regions of bonded 
molecules, at lower T the system is trapped in a metastable excited state where 
large regions of ordered bonded molecules, with different orientations, compete 
in their growth process (left panel on the bottom row in Fig. [3]). These process 
resembles the arrested dynamics of other model systems for glassy dynamics [42|. 
It is characterized by non exponential correlation behavior with characteristic times 
that are several order of magnitudes larger than those at slightly higher T. In our 
case the correlation time largely exceeds the simulation time of 10^ MCS (Full 
Circles in Fig. [2]). 

5. Conclusions 

We have studied the dynamics of the hydrogen bond network of a monolayer of water 
between hydrophobic plates at partial hydration. We adopt a model for water with 
cooperative interactions [35], [271 US]- The model can reproduce the different scenarios 
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proposed for water [32j. We consider the parameters that give rise to a hquid-hquid 
phase transition ending in a critical point [291 l2H]- It has been shown that in the 
supercooled regime, below the TMD line, the system has slow dynamics [301 El l3T] . 
Here we adopt a diffusive MC dynamics to study the case at partial hydration. 

Our results show a complex phenomenology. At high P the effect of the hydrogen 
bond network is negligible due to the high increase of enthalpy for the hydrogen bond 
formation. This is reflected in the dynamics by a moderate increase of the correlation 
time also at very low T. Therefore, the system reaches the equilibrium easily. As a 
consequence, we typically observe the formation of a single large dehydrated region. 

At a pressure close to the liquid-liquid critical pressure, the rapid formation of the 
hydrogen bond network determines the presence of large amount of heterogeneity in 
the system. As a consequence the correlation function shows a strong non-exponential 
behavior at low T, where we typically observe a well developed hydrogen bond network 
with many local heterogeneities. Nevertheless, the hydrogen bond network is fully 
developed only at very low T, keeping the correlation time short and allowing the 
formation of large dehydrated regions on the hydrophobic surface. 

At lower P, the number of hydrogen bonds grows in a gradual and more 
homogeneous way when the temperature decreases. This process slows down the 
correlation decay of orders of magnitudes. Nevertheless, the effect on the exponential 
character of the correlation function is small, because the system has only a small 
amount of heterogeneity. However, the effect on the process of dehydration of the 
surface is evident. 

This effect is extremely strong at very low P, where regions of ordered bonded 
molecules with different orientation grow in a competing fashion. This process give 
rise to an arrested dynamics, typical of many glassy systems, with the correlation time 
that increases of many orders of magnitude for a small decrease in temperature. The 
dehydration process in this case is almost completely frozen on the time scale of our 
simulations. 
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Figure 1. Density p of fluid water as function of temperature T at different values of 
pressure P. (a) At high T a discontinuity in p appears below Pv^/e = 0.5, marking 
the first order phase transition between gas (at low p) and liquid (at high p), ending 
in a critical point at approximately ksTc/f- = 1-9 ± 0.1, Pcva/e — 0.50 ± 0.05. In 
the liquid phase, isobars show a temperature of maximum density (TMD line, marked 
by the thick black line) one of the characteristic features of water, (b) At low T the 
density rapidly decreases for P < J/vhb = le/wo with a discontinuity marking the 
first order phase transition between high density liquid (at high p) and low density 
liquid (at low p), ending in a critical point at approximately ksTc'/c — 0.06 ± 0.01, 
Pc'Vo/e = 0.60 ±0.15. In both panels, symbols are results of the MC calculation, lines 
are guides for the eyes, errors are comparable to the size of the symbols, labels on the 
right of each panel are the values of Pvo/e of the nearest curve. 
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Figure 2. The correlation function CM{t) for different pressure and temperature below 
the temperature of maximum density (TMD). The correlation decays as an exponential 
(dotted lines) at any T > OAe/ks, just below the TMD, for any P < l.Qe/vQ (green 
symbols) and for all T aX P — le/wg (Triangles). The correlation CM{t) is well 
described by a stretched exponential (see text) for P = 0.7e/vo (Diamonds) and 
P = OAe/vo (Squares) for T < O.le/fcs with a stretched exponent f3 that is larger 
at lower T. At low pressure P = Q.le/vo (Circles), CM{t) is non-exponential only at 
very low T = O.OSe/fcs (fuU Circles). 
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Figure 3. Hydrogen bonds configurations sliowing the liydrogen bond ordering at low 
T and at P < J/vhb = ^(/vhb- Symbols: Turquoise dots represent cells containing 
water molecules with at least one bond free; colored lines represent hydrogen bonds 
between molecules in n. n. cells, with six possible colors corresponding to the (7 = 6 
possible values of the orientational state. Every row corresponds to a different P: from 
top, P = le/wo, P = Q.7e/vQ, P = QAe/vq, P = Q.le/vQ. Every column corresponds to 
a different T: from right, T = 0.4e/fcs, T = O.le/ZcB, T = O.OSe/fcs. The size of the 
box is not in scale with the total volume of the system at the different values of P and 
T. 



